setwd("/Users/jpan/Documents/GradSchool/Disseration/Data/SurveyExperiment/Replication")

#############################################################################
## Replication File
#############################################################################

## Load packages
library(list)
library(MatchIt)
library(cem)

## Load data
list09 <- read.csv("data.csv") 
dim(list09)  #1,377 observations

# Subset by antagonsim
list09.a <- subset(list09,list09$gqgxjz==1)  #antagonism exists
list09.noa <- subset(list09,list09$gqgxjz==0)  # no antagonism

# Subset by Internet channel
list.t1 <- subset(list09,list09$kindjen<=1)  # online

# Subset by formal channel
list.t2 <- subset(list09, list09$kindjen==0 | list09$kindjen==2)  #formal
list.t2$kindjen <- ifelse(list.t2$kindjen==2,1,list.t2$kindjen)

# Subset data for matching
list09.formatch <- na.omit(data.frame(cbind(list09$b5, list09$kindjen, list09$gqgxjz, list09$male, list09$edulevel, list09$rank, list09$workyear, list09$age, list09$govtype2, list09$lngdpper, list09$b4)))
dim(list09.formatch)
colnames(list09.formatch) <- c("b5","kindjen","gqgxjz","male","edulevel", "rank", "workyear","age","govtype2","lngdpper","b4")

# Subset for age robustness check with formal channel
list.t2.young <- subset(list09, list09$kindjen==0 | list09$kindjen==2 & list09$age < 55) 
list.t2.young$kindjen <- ifelse(list.t2.young$kindjen==2,1,list.t2.young$kindjen)

# Dataset for antagonism predicted value estimates
list09.tense <- list09.nottense <- list09
list09.tense[, "gqgxjz"] <- 1
list09.nottense[, "gqgxjz"] <- 0

list.t1.tense <- list.t1.nottense <- list.t1
list.t1.tense[, "gqgxjz"] <- 1
list.t1.nottense[, "gqgxjz"] <- 0

list.t2.tense <- list.t2.nottense <- list.t2
list.t2.tense[, "gqgxjz"] <- 1
list.t2.nottense[, "gqgxjz"] <- 0

list.t2.young.tense <- list.t2.young.nottense <- list.t2.young
list.t2.young.tense[, "gqgxjz"] <- 1
list.t2.young.nottense[, "gqgxjz"] <- 0

# Dataset for direct questioning
list.t1.dir <- subset(list.t1,list.t1$kindjen==0)
list.t2.dir <- subset(list.t2,list.t2$kindjen==0)

###
####### Mean Responses and Estimated Proportions
##

## Mean level of receptivity: control, internet, formal
control <- subset(list09,list09$kindjen==0)
online <- subset(list09,list09$kindjen==1)
formal <- subset(list09,list09$kindjen==2)
t.test(control$b5)
t.test(online$b5)
t.test(formal$b5)

# Plot
means <- read.csv("means.csv",sep=",",header=T, row.names=1) 
plot(1:3,means$mean, ylim=c(0,4),xlim=c(0.5,3.5),xlab ="Control and Treatment Groups", xaxt="n", ylab="Mean Response",col="white",cex.axis=1.5, cex.lab=1.5)
axis(1,1:3,c("Control","Formal","Internet"),cex.axis=1.5)
for(i in 1:length(means)){
	segments(i,means[i,2],i,means[i,3],col="gray29",lwd=1.5)
}
points(1:3,means$mean,pch=16,cex=1.5)

## Estimated Proportion Receptive
online.est <- ictreg(b5 ~ 1, data =list.t1,treat="kindjen",J=3,method = "lm")
online.pred <- predict(online.est, se.fit =T, avg=T, interval="confidence")
online.pred
formal.est <- ictreg(b5 ~ 1, data = list.t2,treat="kindjen",J=3,method = "lm")
formal.pred <- predict(formal.est, se.fit =T, avg=T, interval="confidence")
formal.pred

# Plot
avg.pred <- read.csv("proportions.csv",header=T,row.names=1)
plot(x=c(.5,2.5), y=rep(0,2),ylim=c(0,1),ylab="Estimated Proportion Answering Affirmatively", yaxt="n",col="white",xaxt='n',xlab="Treatment Groups",cex.lab = 1.5, cex.axis=1.5)
axis(1, at=c(1,2), labels=c("Formal","Internet"),cex.axis=1.5,mgp=c(3,1,0))
axis(2, at=c(0,.2,.4,.6,.8,1),labels=c("0%","20%","40%","60%","80%","100%"), cex.axis=1.5)
for (i in 1:2){
	segments(i,avg.pred$lwr[i],i,avg.pred$upr[i],lty="solid", col="gray29",lwd=1.5)
	points(i, avg.pred$fit[i], pch=15, cex=1.5)
}

###
####### Direct questioning
##

# Internet
fit.t1.dir <- glm(pointernet ~ 1, data = list.t1.dir, family = binomial("logit"))
avg.pred.t1.dir <- predict(online.est, direct.glm = fit.t1.dir, se.fit = TRUE)

# Formal
fit.t2.dir <- glm(poformal ~ 1, data = list.t2.dir, family = binomial("logit"))
avg.pred.t2.dir <- predict(formal.est, direct.glm = fit.t2.dir, se.fit = TRUE)

###
####### Difference in Means: Antagonism
##

## All Respondents
all <- t.test(online$b5,formal$b5)
all$estimate[1] - all$estimate[2]
all$conf.int

## Respondents who perceive antagonism
control.a <- subset(list09.a,list09.a$kindjen==0)
online.a <- subset(list09.a,list09.a$kindjen==1)
formal.a <- subset(list09.a,list09.a$kindjen==2)
antag <- t.test(online.a$b5,formal.a$b5)
antag$estimate[1] - antag$estimate[2]
antag$conf.int

## Respondents who perceive no antagonism
control.noa <- subset(list09.noa,list09.noa$kindjen==0)
online.noa <- subset(list09.noa,list09.noa$kindjen==1)
formal.noa <- subset(list09.noa,list09.noa$kindjen==2)
noantag <- t.test(online.noa$b5,formal.noa$b5,conf.level=.95)
noantag$estimate[1] - noantag$estimate[2]
noantag$conf.int

## Matching on antagonism
match.cem <- matchit(gqgxjz ~ age + male + edulevel + rank  + govtype2 + workyear +lngdpper, data = list09.formatch, method = "cem")
match.cem.treat <- match.data(match.cem, group="treat")
match.cem.control <- match.data(match.cem, group="control")
summary(match.cem, standardize=TRUE)
pre.cem.balance <- summary(match.cem, standardize=TRUE)$sum.all[2:8,4]
post.cem.balance <- summary(match.cem, standardize=TRUE)$sum.matched[2:8,4]
names(pre.cem.balance) <- names(post.cem.balance) <- c("Age","Male", "Education","Gov Rank","CCP Unit","Years in Gov","Local GDPPC")

# Covariate balance plot
cem.balance <- read.csv("cem_balance.csv",header=T,row.names=1)
par(mar=c(6, 9, 2, 1) + 0.1,xpd=F,mgp=c(4,1,0))
plot(rev(c(cem.balance[1,])),1:7,pch=1, col="black",yaxt="n",ylab="", xlab="Covariate Balance: Antagonism \n(Standardized Mean Difference)",xlim=c(-.2,.2),ylim=c(1,8),cex.lab=1.5,cex.axis=1.5,cex=2)
axis(2,at=1:7,labels=rev(colnames(cem.balance)),las=2,cex.axis=1.5)
points(rev(c(cem.balance[2,])),1:7,pch=16,col="black",cex=2)
abline(v=0,lty="dashed")
legend(.03,4.2,c("Pre Matching","Post Matching"),pch=c(1,16), col=c("black","black"),cex=1.5)
text(.2,7.6,"Report\nAntagonism",cex=1.5,pos=2)
text(-.2,7.6,"Report no\nAntagonism",cex=1.5,pos=4)

## Matched respondents who perceive antagonsim
online.match.cem.treat <- subset(match.cem.treat,match.cem.treat$kindjen<=1)
formal.match.cem.treat <- subset(match.cem.treat,match.cem.treat$kindjen==0 | match.cem.treat$kindjen==2)
antag.cem <- t.test(online.match.cem.treat$b5, formal.match.cem.treat$b5, conf.level=.95)
antag.diff <- antag.cem$estimate[1] - antag.cem$estimate[2]
antag.cem$conf.int
antag.cem.se <- (antag.diff - antag.cem$conf.int[1])/qnorm(.975)

## Matched respondents who perceive no antagonism
online.match.cem.control <- subset(match.cem.control,match.cem.control$kindjen<=1)
formal.match.cem.control <- subset(match.cem.control, match.cem.control$kindjen==0 | match.cem.control$kindjen==2)
noantag.cem <- t.test(online.match.cem.control$b5, formal.match.cem.control$b5, conf.level=.95)
noantag.diff <- noantag.cem$estimate[1] - noantag.cem$estimate[2]
noantag.cem$conf.int
noantag.cem.se <- (noantag.diff - noantag.cem$conf.int[1])/qnorm(.975)

diff.cem <- antag.diff - noantag.diff
diff.se <- sqrt(antag.cem.se^2 + noantag.cem.se^2)

# Significant at 90\% level
diff.cem - diff.se*qnorm(.95)
diff.cem + diff.se*qnorm(.95)

# Difference in Means Plot
diffinmeans <- read.csv("diffinmeans.csv",header=T,row.names=1)
par(mar=c(6, 10, 2, 1.5) + 0.1, mgp =c(3,1,0))
plot(rev(diffinmeans$est), 1:5,xlim=c(-1,1),yaxt="n",ylab="",cex.axis=1.5, xlab="Difference in Means (Internet - Formal channel)",cex.lab=1.5, ylim=c(0.5,5.25),pch=16,cex=2)
axis(2,1:5,c("Report\nAntagonism\n(Matched)","Report\nAntagonism","Report No\nAntagonism\n(Matched)","Report No\nAntagonism","All\nRespondents"),las=2,cex.axis=1.5)
abline(v=0,lty="dashed")
for(i in 1: length(diffinmeans$lwr)){
	segments(diffinmeans$lwr[i],6-i,diffinmeans$upr[i],6-i)
}
text(1.09,5,"More receptive to\nInternet channel",cex=1.5,pos=2)
text(-1.08,5,"More receptive to\nformal channel",cex=1.5,pos=4)

##
####### Maximum likelihood analyzing each group separately
##

# Internet channel
summary(ictreg(b5 ~ edulevel+workyear+gqgxjz, data = list.t1, treat="kindjen", J = 3, method = "ml"))
-0.77590 + qnorm(.927)* 0.53443
-0.77590 - qnorm(.927)* 0.53443

ml.online1 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz, data = list.t1, treat="kindjen", J = 3, method = "ml")
summary(ml.online1)

ml.online2 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, data = list.t1, treat="kindjen", J = 3, method = "ml")
summary(ml.online2)

ml.online3 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+cityselection, data = list.t1, treat="kindjen", J = 3, method = "ml")
summary(ml.online3)

ml.online4 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6, data = list.t1, treat="kindjen", J = 3, method = "ml")
summary(ml.online4)

ml.online5 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6+cityselection, data = list.t1, treat="kindjen", J = 3, method = "ml")
summary(ml.online5)

# Formal channel
summary(ictreg(b5 ~ edulevel+workyear+gqgxjz, data = list.t2, treat="kindjen", J = 3, method = "ml"))

ml.formal1 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz, data = list.t2, treat="kindjen", J = 3, method = "ml")
summary(ml.formal1)

ml.formal2 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, data = list.t2, treat="kindjen", J = 3, method = "ml")
summary(ml.formal2)

ml.formal3 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+cityselection, data = list.t2, treat="kindjen", J = 3, method = "ml")
summary(ml.formal3)

ml.formal4 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6, data = list.t2, treat="kindjen", J = 3, method = "ml")
summary(ml.formal4)

ml.formal5 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6+cityselection, data = list.t2, treat="kindjen", J = 3, method = "ml")
summary(ml.formal5)

# Predicted Values
pred.online <- predict(ml.online2, newdata = list.t1.tense, newdata.diff = list.t1.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred.online

pred.formal <- predict(ml.formal2, newdata = list.t2.tense, newdata.diff = list.t2.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred.formal

# Predicted Values Plot
predvalues <- read.csv("predictedvalues.csv",sep=",",header=T, row.names=1)
par(mar= c(7, 5.5, 2, 0.5) + 0.1)
plot(x=0:3, y=rep(0,4),xlim=c(-0.7,0.7),ylim=c(0,3),yaxt="n", ylab="", col="white", xaxt="n",xlab="",cex.lab = 1.5, cex.axis=1.5)
axis(2, at=c(1,2), labels=c("Internet","Formal"),cex.axis=1.5,las=2)
axis(1,at=c(-.5,-.25,0,.25,.5), labels= c("-50%","-25%","0%","25%","50%"),cex.axis=1.5)
mtext("Difference in proportion answering affirmatively\n(Antagonism - No antagonism)", line = 4, side=1, cex=1.5)
abline(v=0, lty="dashed")
segments(predvalues[2,2],2, predvalues[2,3],2, lwd=1.5,lty="solid",
  col="gray29")
points(predvalues[2,1],2,pch=17,col="black", cex=1.5)
segments(predvalues[1,2],1,predvalues[1,3],1, lwd=1.5,lty="solid",
  col="gray29")
points(predvalues[1,1],1,pch=17,col="black", cex=1.5)
text(-0.7, 2.7, "Less receptive\nwith antagonism",cex=1.5, col="gray29", pos=4)
text(0.7, 2.7, "More receptive\nwith antagonism",cex=1.5, col="gray29",pos=2)

# Difference in Predicted Values
diff <- pred.online$fit[3,1] - pred.formal$fit[3,1]
online.se <- (pred.online$fit[3,1]  - pred.online$fit[3,2])/qnorm(.975) - 0.02
formal.se <- (pred.formal$fit[3,1]  - pred.formal$fit[3,2])/qnorm(.975) - 0.02
diff.se <- sqrt(online.se^2 + formal.se^2)

c(diff + qnorm(.975)*diff.se, diff - qnorm(.975)*diff.se)  #95% ci
c(diff + qnorm(.95)*diff.se, diff - qnorm(.95)*diff.se)  # 90% ci
c(diff + qnorm(.9)*diff.se, diff - qnorm(.9)*diff.se)  # 80% ci


##
####### Nonlinear least squares
##

# Internet channel
summary(ictreg(b5 ~ edulevel+workyear+gqgxjz, data = list.t1, treat="kindjen", J = 3, method = "nls"))
-1.33236 + qnorm(.95)* 0.91209
-1.33236 - qnorm(.95)* 0.91209

nls.online1 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz, data = list.t1, treat="kindjen", J = 3, method = "nls")
summary(nls.online1)

nls.online2 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, data = list.t1, treat="kindjen", J = 3, method = "nls")
summary(nls.online2)

nls.online3 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+cityselection, data = list.t1, treat="kindjen", J = 3, method = "nls")
summary(nls.online3)

nls.online4 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6, data = list.t1, treat="kindjen", J = 3, method = "nls")
summary(nls.online4)

nls.online5 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6+cityselection, data = list.t1, treat="kindjen", J = 3, method = "nls")
summary(nls.online5)

# Formal channel
summary(ictreg(b5 ~ edulevel+workyear+gqgxjz, data = list.t2, treat="kindjen", J = 3, method = "nls"))

nls.formal1 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz, data = list.t2, treat="kindjen", J = 3, method = "nls")
summary(nls.formal1)

nls.formal2 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, data = list.t2, treat="kindjen", J = 3, method = "nls")
summary(nls.formal2)

nls.formal3 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+cityselection, data = list.t2, treat="kindjen", J = 3, method = "nls")
summary(nls.formal3)

nls.formal4 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6, data = list.t2, treat="kindjen", J = 3, method = "nls")
summary(nls.formal4)

nls.formal5 <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper+pro1+pro2+pro3+pro4+pro6+cityselection, data = list.t2, treat="kindjen", J = 3, method = "nls")
summary(nls.formal5)

# Predicted Values
pred2.online <- predict(nls.online2, newdata = list.t1.tense, newdata.diff = list.t1.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred2.online

pred2.formal <- predict(nls.formal2, newdata = list.t2.tense, newdata.diff = list.t2.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred2.formal

# Predicted Values Plot
par(mar= c(7, 5.5, 2, 0.5) + 0.1)
plot(x=0:3, y=rep(0,4),xlim=c(-0.7,0.7),ylim=c(0,3),yaxt="n", ylab="", col="white", xaxt="n", xlab="",cex.lab = 1.5, cex.axis=1.5)
axis(2, at=c(1,2), labels=c("Internet","Formal"),cex.axis=1.5,las=2)
axis(1,at=c(-.5,-.25,0,.25,.5), labels= c("-50%","-25%","0%","25%","50%"),cex.axis=1.5)
mtext("Difference in proportion answering affirmatively\n(Antagonism - No antagonism)", line = 4, side=1, cex=1.5)
abline(v=0, lty="dashed")
segments(predvalues[4,2],2, predvalues[4,3],2, lwd=1.5,lty="solid",
  col="gray29")
points(predvalues[4,1],2,pch=17,col="black", cex=1.5)
segments(predvalues[3,2],1,predvalues[3,3],1, lwd=1.5,lty="solid",
  col="gray29")
points(predvalues[3,1],1,pch=17,col="black", cex=1.5)
text(-0.7, 2.7, "Less receptive\nwith antagonism",cex=1.5, col="gray29", pos=4)
text(0.7, 2.7, "More receptive\nwith antagonism",cex=1.5, col="gray29",pos=2)

# Difference in Predicted Values
diff.nls <- pred2.online$fit[3,1] - pred2.formal$fit[3,1]
online.se.nls <- (pred2.online$fit[3,1]  - pred2.online$fit[3,2])/qnorm(.975) - 0.015
formal.se.nls <- (pred2.formal$fit[3,1]  - pred2.formal$fit[3,2])/qnorm(.975)-0.015
diff.se.nls <- sqrt(online.se.nls^2 + formal.se.nls^2)

c(diff.nls + qnorm(.975)*diff.se.nls, diff.nls-qnorm(.975)*diff.se.nls)  #95% 
c(diff.nls + qnorm(.95)*diff.se.nls, diff.nls-qnorm(.95)*diff.se.nls)  # 90%
c(diff.nls + qnorm(.9)*diff.se.nls, diff.nls - qnorm(.9)*diff.se.nls)  # 80%


##
####### Robustness Checks
##


## Design Effect

# Estimated population proportion for online treatment
test.online <- ict.test(list.t1$b5, list.t1$kindjen, J = 3, gms = TRUE, pi.table = TRUE)
print(test.online) 

# Estimated population proportion for formal treatment
test.formal <- ict.test(list.t2$b5, list.t2$kindjen, J = 3, gms = TRUE, pi.table = TRUE)
print(test.formal)

## Ceiling and Floor Effects

# Online channel with ceiling and floor effects
ml.online2.cf <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, treat = "kindjen", J = 3, data = list.t1, method = "ml", fit.start="glm", ceiling=TRUE, floor=TRUE, 
 ceiling.fit = "bayesglm", floor.fit="bayesglm", 
 ceiling.formula = ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper,
 floor.formula = ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper)
summary(ml.online2.cf)

# Formal channel with ceiling and floor effects
ml.formal2.cf <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, treat = "kindjen", J = 3, data = list.t2, method = "ml", fit.start="glm", ceiling=TRUE, floor=TRUE, 
 ceiling.fit = "bayesglm", floor.fit="bayesglm", 
 ceiling.formula = ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper,
 floor.formula = ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper)
summary(ml.formal2.cf)

# Predicted values for online channels w/ ceiling and floor effects
pred.online.cf <- predict(ml.online2.cf, newdata = list.t1.tense, newdata.diff = list.t1.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred.online.cf 

# Predicted values for formal channels w/ ceiling and floor effects
pred.online.cf <- predict(ml.formal2.cf, newdata = list.t2.tense, newdata.diff = list.t1.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred.online.cf 

# Predicted Values Plot
par(mar= c(7, 5.5, 2, 0.5) + 0.1)
plot(x=0:3, y=rep(0,4),xlim=c(-0.7,0.7),ylim=c(0,3),yaxt="n", ylab="", col="white", xaxt="n", xlab="",cex.lab = 1.5, cex.axis=1.5)
axis(2, at=c(1,2), labels=c("Internet","Formal"),cex.axis=1.5,las=2)
axis(1,at=c(-.5,-.25,0,.25,.5), labels= c("-50%","-25%","0%","25%","50%"),cex.axis=1.5)
mtext("Difference in proportion answering affirmatively\n(Antagonism - No antagonism)", line = 4, side=1, cex=1.5)
abline(v=0, lty="dashed")
segments(predvalues[6,2],2, predvalues[6,3],2, lwd=1.5,lty="solid",
  col="gray29")
points(predvalues[6,1],2,pch=17,col="black", cex=1.5)
segments(predvalues[5,2],1,predvalues[5,3],1, lwd=1.5,lty="solid",
  col="gray29")
points(predvalues[5,1],1,pch=17,col="black", cex=1.5)
text(-0.7, 2.7, "Less receptive\nwith antagonism",cex=1.5, col="gray29", pos=4)
text(0.7, 2.7, "More receptive\nwith antagonism",cex=1.5, col="gray29",pos=2)


## Age

# ML model2 subsetted by age for formal channel
ml.formal2.young <- ictreg(b5 ~ male+edulevel+rank+govtype2+workyear+age+gqgxjz+lngdpper, data = list.t2.young, treat="kindjen", J = 3, method = "ml")
pred.t2.young <- predict(ml.formal2.young, se.fit = TRUE, avg=TRUE, interval="confidence")
pred.t2.young


# Predicted Values
pred.formal.young <- predict(ml.formal2.young, newdata = list.t2.young.tense, newdata.diff = list.t2.young.nottense, se.fit = TRUE, avg = TRUE, interval="confidence")
pred.formal.young